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Air plasma radiation in Local Thermodynamic Equilibrium (LTE) within cylindrical 
geometries is studied with an application towards modeling the radiative transfer inside 
arc-constrictors, a central component of constricted-arc arc jets. A detailed database of 
spectral absorption coefficients for LTE air is formulated using the NEQAIR code developed 
at NASA Ames Research Center. The database stores calculated absorption coefficients 
for 1,051,755 wavelengths between 0.04 jim and 200 pm over a wide temperature (500 K to 
15 000 K) and pressure (0.1 atm to 10.0 atm) range. The multi-group method for spectral re- 
duction is studied by generating a range of reductions including pure binning and banding 
reductions from the detailed absorption coefficient database. The accuracy of each reduc- 
tion is compared to line-by-line calculations for cylindrical temperature profiles resembling 
typical profiles found in arc-constrictors. It is found that a reduction of only 1000 groups 
is sufficient to accurately model the LTE air radiation over a large temperature and pres- 
sure range. In addition to the reduction comparison, the cylindrical-slab formulation is 
compared with the finite-volume method for the numerical integration of the radiative flux 
inside cylinders with varying length. It is determined that cylindrical-slabs can be used to 
accurately model most arc-constrictors due to their high length to radius ratios. 

I. Introduction 

T he radiative heat flux inside the constrictor of a constricted-type arc heater can have a substantial effect 
on the overall conditions of the flow field inside the constrictor as well as that exiting the nozzle into the 
test section. It is therefore important that this radiation field be accurately predicted and fully coupled with 
the governing equations of the flow. However, the determination of the radiative heat flux is computationally 
expensive due to the highly non-gray spectral properties of the test gas. In order to fully couple the radiation 
field to the flow field, a reduced radiation model must first be developed in order to dramatically lower the 
computational costs and make the coupling practical. 

To this end, several reduction strategies have been developed over the years. Multi-band or binning 
methods have been developed which group individual spectral opacities into contiguous groups based on 
frequency ranges and define a single averaged opacity and source for each range. These methods can be 
very accurate, however they generally require a significant number of groups (i.e. 0(lO 4 )) to approach 
line-by-line (LBL) accuracy. Though this is still a significant improvement over the LBL method in terms 
of computational cost, it is still impractical to use multi-group methods for 2- or 3-D flow calculations. 

The opacity distribution function (ODF) method (sometimes referred to as opacity binning or multi- 
band method) was also developed as a technique for reducing the dimensionality of the wavelength space. 1 
As is done in the multi-band methods, spectral lines are grouped and the average spectral properties are 
used to evaluate approximate solutions of the radiative transfer equation (RTE). However, unlike wavelength 
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banding, the binning approach groups spectral properties based on contiguous opacity ranges and defines 
the average opacity and source weighted by some function of opacity. Here, individual wavelengths are not 
grouped in contiguous ranges but are divided into each group (or bin in this case) based on their average 
opacity value over the range of conditions for a given line-of-sight. ODF methods can be much more accurate 
(per group) than the simple banding method and tend to converge after 0(10) bins, 2 however as with the 
multi-band method, reasonable accuracy is generally limited to relatively homogeneous areas of the flow 
field. Moreover, the methodology for choosing the exact binning strategy is not well defined and the success 
of this method is highly dependent on how the bins are chosen. 

The Planck-Rosseland-Gray (PRG) model was originally developed by Sakai et al. 3 for radiating shock 
layers and later applied to radiation modeling in cylindrical media such as arc-heated flows. 4,5 The basis for 
the PRG model is that many wavelengths can be separated into either the optically thin or optically thick 
regimes for some characteristic length. For all other wavelengths, the model assumes a gray gas absorption 
coefficient. For optically thin and thick wavelengths the Planck and Rosseland approximations, respectively, 
provide a closed form solution of the RTE. The standard RTE is then solved assuming a gray gas mean 
absorption coefficient for the remaining wavelengths. The PRG model has proved successful for a wide range 
of conditions, however two selection criteria which determine the boundaries between the optically thin and 
thick regimes must be periodically determined iteratively in order to match the radiative flux determined by 
more detailed band models. This is a major drawback of the PRG model because it must be tuned to each 
experimental condition and cannot be used as a predictive tool in arc-constrictor design studies. 

In this study, the multi-band and ODF methods are used to create a variety of reductions to spectral 
radiative properties for equilibrium air plasma. The accuracy of the reductions are compared to line-by- 
line calculations in cylinders with simulated arc-constrictor temperature profiles using the cylindrical-slab 
formulation for radiative transfer. In addition, the cylindrical-slab formulation is compared to the finite- 
volume method to determine its range of validity for arc-constrictor geometries. 

II. Numerical Methods 

The flow inside the arc column of a full-scale constricted arc-heater is assumed to be in chemical equi- 
librium due to the relatively high pressures and temperatures and long residence times. For this study we 
will further assume that the flow exists in a thermal equilibrium state as well in order to make tabulation 
easier, though all methods employed could easily be extended to include thermal nonequilibrium effects. 
Under these assumptions, any spectral property of the gas, £a, is uniquely determined by some function of 
the mixture temperature and pressure given the element fractions of the mixture. 

Zx=b(T,P) (1) 


The computation of radiative heat flux inside non-gray participating media is covered in numerous texts 6, 7 
and is only summarized here. For a gas in local thermodynamic equilibrium (LTE), the spectral intensity 
I\ along a ray s is governed by the radiative transfer equation neglecting scattering 

dh 


ds 


— ^A (Ib\ — I\) j 


( 2 ) 


where k\ is called the spectral absorption coefficient and Ib\ is the blackbody emission function or Planck 
function 

1 ^= 2 f h r: . ( 3 ) 

ex P (xkr) ~ 1 

The total radiative heat flux q at a given location is then given by integrating the spectral intensity over 
all wavelengths and solid angles, 


q{x) = 


I\{x , s)dtid\. 


( 4 ) 


A. Absorption Database 

The spectral absorption coefficient, is determined by considering the various radiative processes occurring 
in a high temperature gas; these are typically classified as bound-bound, bound-free, and free-free processes. 
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Chauveau et al. have recently studied these various phenomena in depth and have provided an excellent 
overview of equilibrium air plasma radiation in Ref. 8. For this study, the Nonequilibrium Air radiation code 
(NEQAIR 9 ) developed at NASA Ames Research Center is used to tabulate spectral absorption coefficients 
over a range of temperatures between 500 K and 15 000 K and pressures between 0.1 atm and 10.0 atm for 
LTE air. NEQAIR contains an extensive line database for Air and C0 2 species. Furthermore, the NEQAIR 
code has been used extensively by other researchers and validated against numerous experimental datasets. 

The species number densities for LTE air, which are required to compute the spectral properties in 
NEQAIR, are computed using the recent curve fits of D’Angola et al. 10 The composition of air is assumed 
to be 80% N 2 and 20% 0 2 by volume. The curve fits of D’Angola et al. were computed using the following 
19 species: 


N 2 , N+ n, n+, n 2+ , n 3+ , n 4+ , 0 2 , 0+ 0 2 , 

In addition to the continuous processes (bound-free 
to generate the (bound-bound) spectra in NEQAIR 


O, O , 0+, 0 2+ , 0 3+ , 0 4+ , NO, NO + , e- 

and free- free), the following band systems were used 


Atomic Systems: 
N O 


Diatomic Systems: 

Nj 1 st neg. 

N 2 1 st and 2 nd pos., BHl and BH2, LBH, Carroll- Yoshino, Worley- Jenkins 
0 2 Schumann-Runge 
NO infrared, 7, /?, S, e, 7’, /3’ 


1. Spectral Range 

The total radiative heat flux is determined via Eq. 4 where the limits of integration over A range over the 
entire set of positive real values. In practice however, the spectral range is finite. For this study, the wave- 
length range was chosen in order to capture all important contributions of emission and absorption 

Wien’s displacement law provides the wavelength where the blackbody emission is maximized for a 
given temperature. 

arg max Rx = 2898 um K (5) 

n\T 

This relationship is shown in Fig. 1 which serves as an indicator of the relative importance of each wavelength 
on the total emission (and absorption) for given temperatures. For this study, a wavelength range of 400 A 
to 200 000A (77 = 1/A = 500 cm -1 to 250 000 cm -1 ) was chosen based on Fig. 1. This choice can be easily 
supported by comparing the numerically integrated Planck function over this range with the analytical value 
of the integral of the Planck function from 0 to 00 

roo 

/ I bX d\ = -T\ (6) 

Jo 7T 

where a is called the Stefan-Boltzmann constant and is defined to be 

97T 5 fr 3 4 

o- = , b = 5.6697 x 10 -5 erg/cm 2 K 4 s. (7) 

15/i d cr 

Fig. 2 shows such a comparison. It is clear from the figure that the wavelength range of 0.04 pm to 20 pm is 
sufficient to capture the dominant emission and absorption phenomena over all temperatures and pressures 
of interest. 

2. Grid Spacing 

The grid spacing used to tabulate the absorption coefficient database was determined through a balance 
between minimizing the grid size while ensuring that the spectral properties are fully resolved over the entire 
wavelength domain. To this end, a simple grid resolution study was performed to determine the optimal grid 
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Wien's Displacement Law 


Validation of Planck Function Integral (10 bins) 




sicjma/pn 
P = lOati 


Temperature [K] 


Figure 1: Wien’s distribution showing affective 

wavelength range for temperatures between 300 K 


and 15 000 K. 


Figure 2: Comparison of the integrated Planck func- 
tion with limits of integration from 0.04 pm to 20 pm 
and the theoretical curve for 0 to oo. 


spacing necessary to meet both criteria. It was determined that the following grid ranges and spacings in 
Tabic 1 provide fully resolved k\ while minimizing erroneous wavelength points. The total size of the entire 
grid is 1,051,755 wavelengths. 


Table 1: Wavelengthgth grid used in neqair.inp. 


Ai [A] 

A 2 [A] 

AA [A] 

Points 

400 

2,000 

0.04 

40001 

2,000 

6,350 

0.10 

43501 

6,350 

200,000 

0.20 

968253 


The resulting absorption coefficient spectrums for five different temperatures at 1 atm for LTE air as 
computed by NEQAIR are presented in Fig. 3. From the figure, it is clear that the spectral absorption 
coefficient varies wildly with wavelength and temperature due to the highly non gray nature of the absorb- 
ing/emitting gas and the changing chemical composition due to the LTE air formulation. Solving the RTE 
in Eq. 2 requires integrating over each spectral wavelength, line of sight, and solid angle. Performing the 
full integration in this way is called a line-by-line calculation and represents the most accurate method for 
computing the radiative transfer in a non gray medium. For 3-dimensional geometries, this can pose an 
impractical computational task, especially when it is desired to couple the radiation field to the flow field 
governed by the Navier-Stokes equations. (Note: the strong, isolated downward-pointing spikes in Fig. 3 
(a), (b), and (c) are artifacts of the wavelength domain decomposition and are not physical.) 


B. Multi-band Model and Opacity Distribution Functions 


In essence, the multi-band model and the opacity distribution function model (multi-bin) are similar methods 
in that they both attempt to reduce the dimensionality of the wavelength space by grouping individual 
wavelength’s together and defining an average absorption coefficient (opacity) and emissivity for each group. 
The primary difference lies in how the groups are formed. This common ground allows a succinct development 
of both methods starting with an identical formulation but diverging in the formation of the groups. 

The radiative heat flux is computed using Eq. 4. Since the direction-cosine of the ray s is not a function 
in wavelength (or any spectral quantity), the integral over wavelength in Eq. 4 may be moved inside the 
radiative transfer equation, Eq. 2. 


t 3 

ds 


I x d\ = / KxhxdX - 


K X I\dX 


( 8 ) 
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Absorption Coefficient [cm' 1 ] Absorption Coefficient [cm’ 1 ] 



25000 50000 75000 100000 125000 150000 

Wavenumber [cm’ 1 ] 



25000 50000 75000 100000 125000 150000 

Wavenumber [cm' 1 ] 


(a) 5000 K 


(b) 8000 K 



(c) 10 000 K 


(d) 12 000 K 



(e) 15 000 K 

Figure 3: Absorption coefficient spectra for an LTE air plasma at 1 atm and multiple temperatures. 


Note that the above equation does not make any assumptions or approximations and is still completely valid 
in a mathematical sense. However, the right-hand term of the right side of Eq. 8 is not closed and cannot 
be determined. In order to close the equation the following approximation is made. 

/o°° n x I x dX ~ f 0 °° Kx l bx d\ 

fo°° I x dX ~ r h\dX 1 j 
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from which it follows that 



K\I\dX 


The integrated RTE can then be written as 


Jo K\h\d\ 

r hxdx 





(10) 


( 11 ) 


where 

r°° r°° f n °° K\i\dx , s 

1= I\dX, I h = I I bx dX , and re = J ° * . ( 12 ) 

Jo Jo 1 b 

In its present form, Eq. 11 can be solved using any RTE solver, and the corresponding heat flux is then 



d I[x) dtt 


(13) 


Furthermore, since n and I b can be tabulated for a range of temperatures and pressures, the integrated 
RTE is extremely cheap to solve - reducing O(l0 6 ) equations to 1. However, as the absorption coefficient is 
rapidly varying function of wavelength, the Planck-weighted absorption coefficient k is often very inaccurate 
in representing the intensity- weighted function (the approximation made in Eq. 10). In order to improve this 
accuracy, the integrated RTE is solved over subsets of the entire spectral range using multiple Planck-mean 
absorption coefficients which are either spectrally local (banding) or local to a given opacity range (binning) . 

The multi-band models fall into the first category of spectrally-local mean absorption coefficients. In 
these models, the wavelength range is divided into contiguous sets of wavelengths called bands. 


{Ai} MB = {A € R + | Xf < X < Af } V i = 1 . . . riband 


(14) 


A natural band sizing comes from a logarithmic scaling as spectral features tend to be much finer at low 
wavelengths (absorption/emission due to changing electronic levels) and more spread out at large wavelengths 
(continuous and rovibrational band radiation). 

So-called opacity distribution functions are found if the grouping is based instead on contiguous ranges 
in absorption coefficient space, where each range is called a bin. The wavelengths included in each bin i are 
then described by 

(AJodf = {A e R+ I K i ^ K (A) < ^i 7 } V i = 1 . . . n bin (15) 

Note that this grouping rule will not generally produce contiguous wavelength segments. As with the bands, 
the bin levels are scaled logarithmically. Once the actual wavelength groupings are determined the integrated 
RTE is written for each group (bin or band) 


9Ii_ . (T 

\lbi J-i) 5 

where the integrals in Eq. 12 are evaluated over each bin or band 


h = 


1 


'{Ai} 


I\d\, I b i = / h\dX, and k* = — / K\I b \dX, 


'OH 


'{Ad 


(16) 


(17) 


and the total intensity is simply the summation of the intensities from each group, I(x) = ]Tb Ii(x). 

The multi-group method 1 combines the multi-bin and multi-band methods into a general framework. This 
is done by considering first a set of contiguous bands which split the spectral range into smaller segments 
which are then further split into separate bins leading to n groU p = ^band^bin total groups over which the 
absorption coefficient is averaged. For this study, a multi-group approach is taken. Note that the multi-band 
and -bin methods are simply special cases of the multi-group method where either only one bin per band or 
one band total are used. 

The procedure for generating the groups from the tabulated spectral absorption coefficients is as follows. 
First, the spectral range is split into n ban< i number of bands using a logarithmic scale in wavelength. This is 
done in order to try and even out the number of wavelengths associated with each band as the wavelength 
grid spacing increases with increasing wavelength. Another important reason is that a primary assumption 
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400 

Wavelength [nm] 


(a) Unordered absorption coefficient spectrum. 



(b) Ordered absorption coefficient spectrum. 


Figure 4: Example showing how the binning procedure works with 4 bins in the 300 nm to 500 nm range for 
equilibrium Air at 10 000 K and 1 atm. 


of the binned RTE methodology is that the Planck function is considered constant over each group. From 
Fig. 1 it is clear that Ib\ varies exponentially at lower wavelengths. In addition to the logarithmic scaling, 
an additional constraint is placed on the bands. It was found during the course of testing that since above 
100 pm, n\ is a smooth function of A, it was beneficial to force the last band to cover the 100 pm to 200 pm 
range. This constraint improves accuracy by allowing the remaining bands to cover a smaller spectral range. 

After the wavelengths are split into bands, each band is then subdivided into ribin bins. Several strategies 
for determining the opacity ranges which determine the bins in each band were tested. If a constant or 
logarithmic spacing was used, it is was found that not all bins were guaranteed to contain at least a single 
wavelength. In addition, it was likely that one bin would contain over 90% of the wavelengths in the band 
and thus the other bins would suffer from a statistical point of view. In light of these issues, another method 
was adopted in which the spectral absorption coefficients were first ordered in each band and the bin ranges 
were determined by assigning equal numbers of wavelengths to each bin. This process is shown graphically 
in Fig. 4. 

After assigning each wavelength to a group, the mean absorption coefficient and integrated Planck func- 
tion as defined in Eq. 17 were computed for each group and tabulated versus temperature and pressure to 
form a reduced spectral model. One important observation about the integration method should be noted. 
Originally, the authors attempted to use the trapezoidal rule for computing the integrals in Eq. 17, however, 
as can be easily verified, the trapezoidal rule is not consistent with the original RTE. In other words, when 
the number of groups approaches the number of wavelengths, the accuracy of the “reduced” model should 
converge to LBL, however this is not the case when trapezoidal rule is used because cross-multiplication 
terms are introduced between spectral values. These inconsistencies actually increase the error of a reduc- 
tion substantially as the number of groups increase. Therefore, the rectangular rule was used instead to 
perform the numerical integration which is consistent with LBL calculations. 

C. Cylindrical Slab Formulation 

The equations for radiative heat flux in a cylindrical “slab” are derived in the paper by Kesten et. al. 11 To 
compute the radiation intensity in the radial direction at a point O (see Fig. 5), they analytically integrate 
the equation representing the change in intensity along a light beam, Eq. 2, starting at a point A, on the 
cylinder wall, passing through point O. 

After a series of coordinate transformation and analytical integration over angles, they write, 


«(r) =4 /; /2 [/,UD 3 + jV*--*. y KnuW 


w 


r-sj R-r 2 sin 2 7 „ 

+ Jo 


{ /o’ COS7 Knu(y")dy" + K n u{y")dy"} r]i(y')dy' 
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Figure 5: Ray path (AOB) in cylindrical geometry. The coordinates (s and x) and angles are shown on the 
figure. 
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Figure 6: Ray paths during integration in 7. Two rays are shown, one corresponding to A’OB’ (0 < 7' < 7r/2), 
the other corresponding to A” OB” with 7" = 7r — 7'. 


+jr' lD 2{j^ 7 <s"W} n,WW 

~h\ rd 3 


cos 7^7 


K ' 2 * 


+/, 


n,(y"W‘ 


r cos 7 


cos 7^7 


(18) 


where Ii\r is the radiative intensity at the wall and D n is the exponential integral function, 

r»l ,,n— 1 


D n {x) = / ; = exp {-xl\i)d[i 

Jo — 


Nicolet et al. 12 developed a numerical algorithm to solve Eq. 18. They start by adopting an approxima- 
tion to D n using exponential fits of the form, 13 

D n (x) = a n exp(— & n x), (19) 

and substitute the above in Eq. 18. They chose a best fit to D 3 , and have set, 

b 2 = b 3 = 5/4 

a 2 = a 3 = tt/4 

The cylinder is then discretized into N mesh points. They derived the final approximation to Eq. 18 
using logarithmic quadrature rules to obtain, 

%i = - %i) 
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where, 


9m = 4 E 

J=2 


AT 


G 


i,j — 1 


^7-1 


(20) 


Nicolet et al. 12 used a band approximation where they have set the radiative energy flux at the 
wall, to be the integral of the emissivity in a band [h,h+ 1 ], determined by, 


h\R = <TT*(m)-m + 1 )) 


where /(£;) is the fractional function, 



g 3 

exp(0 - 1 


dd 


and ^ = hc/(Xik B T). 

The expressions for GT- and Gfj in Eq. 20, are Eq. (A-21) and (A-22), given in Nicolet et al. 12 

D. Finite Volume Radiative Transfer 

For verification of the cylindrical slab results, the radiative transfer problem was also solved by means of a 
3-d finite-volume method using the same spectral absorption coefficients and emission function. We write 
the radiative transfer equation, Eq. 2, in three-dimensional form as 


■ V/ A (x, 17) = /t A (x) (I b a(x) - I A (x, 17)) , (21) 

where 1 l is the direction of propagation and x is the spatial position. The distance variable s in Eq. 2 is 
measured along the vector 17. Upon integration of Eq. 21 over a computational cell we obtain 


' cell surface 


/a(x, 17)17 • dS 


re A (x) (Aa(x) - Ja(x, 17 )) dV. 


/cell volume 


(22) 


Discretizing Eq. 22 by approximating the volume integral with the value of the integrand at the cell center 
times the volume and approximating the surface integral by a summation over the faces of face-centered 
values times the corresponding areas, one obtains 


E A(xfc, 17)17 • AS fc = ka(x c ) (/ 6 a(x c ) - /a(x c , 17)) V, (23) 

faces k 


where x c is the position of the cell center, x*, is the position of the center of face k, V is the cell volume, 
and AS fc is the outward-pointing surface area vector of face k of the cell. Splitting the surface summation 
into incoming and outgoing parts and using cell-center values for the outgoing facial intensities in each cell, 
one gets the standard finite-volume method for radiative transfer: 7 


/a(x c , 17 ) 


ka(x c )/{,a(x c )U + X)fc,f2-As fc <o A( x fc, 17) 1 17 • ASfcl 

(x c )U + ASfc>o ^ 


(24) 


While this classic formula, Eq. 24, is an adequate approximation for many applications, we used a similar 
but more accurate formulation. 14 This method allows for significant variation in /a(x, 17) across a cell. The 
outgoing intensity at each outgoing face is given in this method by 


/a(x 0 , 17 ) = 


Sfc,n-As fc <o d\{^-k, O) 1 17 • AS^I _ 


E 


/c,n-AS fc >0 


17 • AS* 


KA(Xc)rf 


(1 




" )d )Ibx(*c), 


where x„ is the position of an outgoing face, and 


(25) 


d = 


V 


(26) 


Efc,f2-AS fc >0 ^ 

is the average value of the distance across the cell from the incoming to the outgoing faces in direction 17. 
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The value of I\(x 0 , ft) is computed using Eq. 25 for a set of ft vectors, the so-called discrete ordinate 
method. ' This set is chosen to produce an accurate quadrature over the 47 t direction space for computing 
the divergence of the radiative heat flux, V • • The average radiative heat flux in a cell is thereby computed 

via 


1 

V 


X7q x dV 


/ cell volume 


1 

v 

i 

V 
1 

V 


Vx-dS 


/ cell surface 



E I x(x k , ft) ft ■ AS k dfl 

faces k 


E E I\(x k - ft,) ft, ■ AS k AQj 

i faces k 


(27) 

(28) 

(29) 


where the weights Af are given by the chosen quadrature formula. In this work all angular quadratures 
were done with the S4 set of discrete ordinates, which uses 24 directions ft t . 


III. Results and Discussion 


A. Reductions 


High pressure arc-columns typically have a hot core near the upstream region which expands and cools 
downstream. The radial temperature profile inside a column resembles a gaussian temperature distribution 
with a peak temperature in the center and a fixed wall temperature at the radius of the column. The 
temperature distribution of Eq. 30 has been created in order to simulate such profiles at different downstream 
cross-sections. 


T(r) 


(T wa u — T max ) 


exp (— (r/i?) 2 /2cr 2 ) - 1 
exp (— l/2cr 2 ) — 1 


+ T m 


(30) 


The temperature distribution is parameterized based on the maximum and wall temperatures, T max and Twain 
the radius of the column, R, and a stretching factor, cr, which corresponds to the variance of the gaussian 
profile. Three temperature profiles have been chosen to test the performance of the spectral reductions using 
Eq. 30. The values of T max , T wa u, R , and <7 are given in Table 2 and the actual temperature profiles are 
plotted graphically in Fig. 7. 


Table 2: Parameters corresponding to model temperature distribution used to test the accuracy of various 
model reductions (see Eq. 30). 


Name 

T 

- 1 - max 

T wa ll 

R (cm) 

a 

G08 

8,000 

1,000 

1.5 

1000 

G10 

10,000 

1,000 

1.5 

1 

G12 

12,000 

1,000 

1.5 

0.5 


Several reductions were computed using the methodology described in Section B. A variety of band and 
bin combinations were chosen in order to study the overall effectiveness of the methodology in accurately 
reducing the full spectrum of absorption coefficients to a manageable number. Each reduction was compared 
to LBL calculations using the cylindrical-slab method for the three temperature distributions detailed in 
Table 2. The accuracy of each reduction is assessed based on two criteria. The first criteria, cv-q, is based 
on the relative, maximum percent difference of V • q(r) as compared to the LBL calculation along the entire 
radial distance in the column. 


e v<? 


m»x r (7 ■ fo, L (r) - V ■ #■)) x m% 
max r |V • <?lbl v ) I 


(31) 


Note that defining the error in the radiative source term in this way does not lead to large relative errors 
when V • q(r) is nearly zero. Instead, this criterion gives a sense of the maximum error with respect to the 
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14000 


G12 

G10 



Radial Location [cm] 

Figure 7: Gaussian temperature profiles from Table 2. 


overall order of magnitude of the peak heating within the column. The second criteria, e qr wall , is the relative 
error in the radiative heat flux at the wall as compared to LBL. 


e 


Qr, wall 


<7r, LBL (-R) ~~ q r (R) 

q r ,hBh(R) 


x 100% 


(32) 


Table 3 summarizes the 28 different reductions computed and their accuracies as compared to LBL based 
on the above two criteria. The table is ordered by increasing bands and then increasing bins to show the clear 
relationship the number of bands has on the overall accuracy of the reduction. Regardless of the number of 
bins used in the reduction, an error of less than 10% for all temperature profiles cannot be obtained with 
less than 1000 bands. It is also important to note that in general, it appears that the radial heat flux at the 
wall is easier to accurately predict than the divergence of the heat flux within the column, however for all of 
the pure binning cases (1-5), the opposite trend is observed. 

Fig. 8 presents the error criterion of Table 3 in a graphical way making it easier to visualize the overall 
trends. It is clear from the figure that indeed the wall radiative heat flux is better predicted by the reduced 
models as compared to the divergence of the flux in the column. The bulk of the e g all values are centered 
around 10% where as the ev-g values trend closer to 100%. Fig. 8 also highlights another important (though 
expected) result. The colors in the figure represent which temperature profile in Table 2 the errors belong 
to. It is therefore clear to see that in nearly all cases, the lower the peak heating is, the better each reduced 
model can predict the radiative heating within the cylinder. This suggests that a practical RTE solver could 
select from a range of reduced models depending on the core temperature of the arc-column, enabling a 
computational speedup in the downstream regions where the temperatures are more uniform. 

Fig. 9 shows the actual q r and V • q profiles for the LBL and five reductions in Table 3 which have the 
lowest overall errors. The heat flux within the column first peaks to a maximum due to the hot core at the 
center before the absorption in the colder regions towards the wall tends to equilibrate the radiative flux from 
both directions. The radiative flux at the wall is much lower than the peak location, owing to the absorption 
in the wall region. Interestingly, the errors in the peak heating region due to the various reductions shown 
in Fig. 9a have little affect on the final radiative flux reaching the wall, suggesting that the cold region is 
optically thick. 

In addition, it is clear that small differences in the core temperature result in very large differences in 
peak heating as can be seen by comparing the V -q profiles for each temperature distribution in Fig. 9b. This 
is most likely due to the fact that at 10 000 K, atomic nitrogen is at its peak concentration before ionization 
begins for equilibrium air at latm. Above 10 000 K, the concentrations of N + and e~ begin to become 
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Table 3: Summary of the errors, ey.q- (Eq. 31) and e 9rwall (Eq. 32), of each reduction tested on the gaussian 
temperature profiles listed in Table 2. 


G08 G10 G12 


No. 

Type 

^band 

^-bin 

^group 

rr \ 

<1 

^ Qr , wall 

ev-q 

^ Qr , wall 


^ Qr , wall 

1. 

bins 

l 

10 

10 

- 26.0 

- 149.7 

- 84.6 

- 122.1 

- 121.8 

- 179.7 

2 . 

bins 

l 

20 

20 

- 26.0 

- 147.9 

- 82.2 

- 119.5 

- 118.5 

- 174.7 

3 . 

bins 

l 

50 

50 

- 25.9 

- 147.0 

- 79.4 

- 117.2 

- 116.4 

- 171.6 

4 . 

bins 

l 

100 

100 

- 25.9 

- 146.5 

- 78.9 

- 116.7 

- 116.2 

- 170.9 

5 . 

bins 

l 

500 

500 

- 25.9 

- 145.4 

- 76.4 

- 114.8 

- 114.2 

- 167.9 

6 . 

multi-group 

4 

50 

200 

- 9.9 

- 9.4 

- 81.4 

- 8.6 

- 113.5 

- 12.6 

7 . 

multi-group 

5 

5 

25 

- 47.5 

- 41.6 

- 227.6 

- 92.4 

- 170.7 

- 115.5 

8 . 

multi-group 

5 

10 

50 

- 43.2 

- 13.8 

- 207.6 

- 21.4 

- 153.4 

- 26.2 

9 . 

multi-group 

5 

20 

100 

- 38.3 

- 8.5 

- 172.0 

- 8.4 

- 129.7 

- 9.2 

10 . 

multi-group 

8 

25 

200 

- 8.3 

- 4.4 

- 72.8 

- 5.3 

- 94.2 

- 6.0 

11 . 

bands 

10 

1 

10 

- 174.5 

- 102.4 

- 733.7 

- 397.1 

- 475.7 

- 482.5 

12 . 

multi-group 

10 

5 

50 

- 122.2 

- 10.8 

- 322.8 

- 16.1 

- 212.5 

- 25.5 

13 . 

multi-group 

10 

10 

100 

- 77.0 

- 7.9 

- 129.2 

- 9.6 

- 100.0 

- 16.2 

14 . 

multi-group 

10 

20 

200 

- 32.6 

- 7.0 

- 50.1 

- 8.3 

- 57.5 

- 13.0 

15 . 

multi-group 

10 

50 

500 

- 3.9 

- 6.6 

- 32.2 

- 6.8 

- 51.3 

- 9.2 

16 . 

multi-group 

20 

5 

100 

- 44.9 

- 6.0 

- 175.4 

- 8.6 

- 134.4 

- 15.2 

17 . 

multi-group 

20 

10 

200 

- 27.8 

- 4.7 

- 99.2 

- 6.4 

- 85.4 

- 11.0 

18 . 

multi-group 

20 

25 

500 

- 12.2 

- 4.3 

- 41.6 

- 4.2 

- 55.8 

- 6.6 

19 . 

multi-group 

25 

8 

200 

- 23.6 

- 3.6 

- 82.6 

- 6.2 

- 69.2 

- 11.3 

20 . 

multi-group 

25 

20 

500 

- 12.0 

- 3.5 

- 39.0 

- 4.8 

- 47.9 

- 7.8 

21 . 

multi-group 

50 

4 

200 

- 26.7 

- 2.6 

- 79.3 

- 4.7 

- 59.1 

- 8.2 

22 . 

multi-group 

50 

10 

500 

- 12.4 

- 2.9 

- 31.6 

- 5.0 

- 34.9 

- 7.8 

23 . 

bands 

100 

1 

100 

- 50.7 

2.6 

- 161.7 

- 15.0 

- 136.0 

- 23.4 

24 . 

multi-group 

100 

2 

200 

- 26.5 

- 1.2 

- 75.2 

- 9.9 

- 58.2 

- 14.8 

25 . 

multi-group 

100 

5 

500 

- 10.5 

- 2.9 

- 27.4 

- 7.3 

- 27.8 

- 10.1 

26 . 

bands 

500 

1 

500 

- 16.5 

4.6 

- 40.6 

- 7.0 

- 38.4 

- 12.9 

27 . 

bands 

1000 

1 

1000 

- 5.8 

4.8 

- 7.3 

- 0.8 

- 3.9 

- 1.9 

28 . 

bands 

10000 

1 

10000 

4.9 

4.1 

- 2.2 

2.0 

1.1 

3.5 


substantial which increase the radiant power due to bound-free and free-free transitions. Thus, the peak 
heating at 12 000 K is about — 16kW/cm 3 while at 10 000 K and 8000 K the peak heating is only — 2kW/cm 3 
and — 0.5kW/cm 3 respectively. This observation underlines the importance of accurately predicting the gas 
composition in the column which can be highly dependent on the radiative-flow coupling. 

Finally, we consider the poor performance of the pure binning reductions (reductions 1-5 in Table 3). 
Surprisingly, the performance of the reductions is barely affected by the number of bins used. For example, 
the error in the radiative flux at the wall for the G12 temperature profile changes from 179.7% with just 10 
bins to only 167.9% with 500 bins. In contrast, the pure banding yields for the same error criteria 482.5% 
error with only 10 bands but just 12.9% with 500 bands. 

Fig. 10 shows the q r and V • q profiles for each temperature distribution using each of the pure binning 
reductions as compared to the LBL calculations. As hinted at by the errors listed in Table 3, the pure 
binning profiles are nearly identical regardless of the number of bins used. Several authors have previously 
noted that the assumption made in Eq. 16 is that each group of wavelengths (either bins or bands or mixed) 
either spans a spectral region in which the Planck function is nearly constant, or the spectral absorption 
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Figure 8: Results of Table 3 shown graphically. Symbols represent type of reduction (0, pure binning; □ , 
pure banding; A, mixed). Colors represent the temperature distribution (red, G12; green, G10; blue, G08). 




Figure 9: Comparison of LBL and reduced model calculations of the radiative heat flux and its divergence 
for the temperature distributions shown in Fig. 7. Only the five reductions with the smallest errors as given 
in Table 3 are shown. 


coefficient can be treated as randomly distributed such that the errors associated with the Planck-weighted 
mean absorption coefficient are small. If either assumption is not met, then the errors in the reduction will be 
large. Therefore, the poor performance of the pure binning reductions as shown in Fig. 10 can be attributed 
to the fact that the bins may span the entire wavelength spectrum, in which case the Planck function can 
range over many orders of magnitude. This is further substantiated by the observation that using a few 
bands to divide the wavelength spectrum into smaller regions can dramatically increase the accuracy of the 
reductions. For instance, using 500 bins alone (reduction 5) yielded a e g n value of -167.9%, however if the 
total wavelength range is split into 10 subregions with 50 bins each (reduction 15), this error drops to -9.2% 
(note that the error in the divergence term is still over 51% for this case). 

B. Cylindrical-Slab vs. Finite- Volume 

In the following section, a comparison is made between the cylindrical-slab formulation and the finite-volume 
formulations detailed in previous sections. This is done in order to both verify the correctness of the the 
cylindrical-slab formulation as well as to determine when it is valid to use this method for arc-jet geometries. 
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Figure 10: Comparison of LBL and pure binning calculations of the radiative heat flux and its divergence 
for the temperature distributions shown in Fig. 7. 


To begin, Eq. 30 is rewritten in two dimensions by linearly interpolating from a constant wall temperature, 
Twain at z = 0 to a gaussian temperature profile with peak temperature, T max , at z = L. 


T(p, 0 


(1 - C)Twall + C 


(T wa ll — T m ax) 


expo jp, C) - 1 
expa(l, C) — 1 



(33) 


where 



jr, a(p, C) 


-0.5p 2 
(C<7l + (1 - 


(34) 


The benefit of such a parameterization is that it is nondimensionalized by the cylinder’s length and radius, 
thus radial temperature profiles at a particular non dimensional cross-section, remain unchanged when 
the length of the cylinder, L, is changed. As with the stretching factor <7 in Eq. 30, the factors <Ti and <72 in 
Eq. 34 determine how stretched the gaussian profiles are at z = L and z — 0 respectively. For the remainder 
of this section, these values are taken to be <j\ = 0.2 and 02 = 1.0. Fig. 11 shows the progression of radial 
temperature profiles for these stretching factors in the non dimensional coordinates of Eq. 34. 

The finite- volume method was used to compute the radiative flux and its divergence in a cylinder whose 
temperature profile is constant in the angular direction and follows the temperature parameterization of Eq. 
33 in the radial and axial directions. The radius of the cylinder, R , was taken to be 3 cm and four separate 
lengths were chosen corresponding to L = 101?, 20 R, 30 R and 40 R. The computational domain was split into 
300 radial cells, 300 axial cells, and 20 angular cells totaling 1.8 million cells to model the entire domain. In 
order to limit computational time, the 1000 band reduction (number 27 in Table 3) was used to represent 
the spectral properties of the gas. 

Fig. 12 presents the radial heat flux and the divergence as a function of the radial location at the 0.5L 
cross-section of the cylinders (averaged over each angular location) compared to the solution obtained for the 
same temperature profile using the cylindrical-slab formation with the 1000 band reduction. As can be seen 
from the figure, as the L/R. ratio increases, the cylindrical-slab formulation becomes closer to reproducing 
the exact integration for the radiative heat flux in the radial direction. This could be expected but it 
is far from guaranteed as the assumptions made by the cylindrical-slab formulation are 1) the cylinder is 
infinitely long, and 2) the spectral properties vary only in the radial direction. The first assumption is clearly 
approached with an increasing L/R ratio, however the second assumption is less obvious. However, from the 
above results, it is clear that allowing small deviations in the radiative properties in the axial direction is a 
sufficient condition for the applicability of the cylindrical-slab formulation. 

Fig. 12 shows that for the particular temperature parameterization given, the cylindrical-slab formulation 
is accurate for cylindrical geometries with L/R > 30. For comparison, Nicolet et al. 12 collected size 
information of various constricted arc heaters in use in 1975. The L/R values for these constrictors ranged 
from about 14 to over 162 where all but one had a minimum L/R of 30. The Interactive Heating Facility 
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Figure 11: Radial temperature profiles given by Eq. 33 at several cross-sections along the cylinder for 
stretching factors u\ = 0.2 and a 2 = 1-0. 




Figure 12: Comparison of the radial heat flux and the divergence of the heat flux for finite- volume calculations 
performed with varying the length of the cylinder versus a cylindrical slab formulation. 


(IHF) and the Aerodynamic Heating Facility (AHF) at NASA Ames Research Center operate today at L/ R 
ratios of 97.3 and 79.6 respectively. 

In terms of performance, the cylindrical-slab formulation was two orders of magnitude faster than the 
finite-volume solver after taking into account that the cylindrical-slab code must be used at every axial 
cross-section. In addition, because each cross-section is treated as independent, the cylindrical-slab lends 
itself to a higher degree of parallelization. This increase in performance along with a complete reduction in 
spectral properties makes fully coupled flow/radiation calculations possible for arc-constrictor simulations 
while maintaining a high level of accuracy. 
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IV. Concluding Remarks 


Full LBL calculations of radiative transfer for real problems is in general computationally impractical 
with todays computational facilities. This is further compounded when coupling the radiative source terms 
with the energy conservation equation for fluid dynamics simulations. For cylindrical geometries, the above 
results show that the cylindrical-slab formulation can reach the same accuracy as performing the full spatial 
integration when the cylinder is long enough. One important application of this methodology is modeling 
the radiative transfer occurring in high pressure arc-constrictors which greatly impacts the flow inside the 
constrictor. The cylindrical-slab formulation provides a considerable increase in performance over the finite- 
volume method while retaining reasonable accuracy. 

The computational performance can be further improved when a reduced radiation model is used. It was 
shown in the above that a 1000 band reduction is sufficient to completely model equilibrium air radiation over 
a wide temperature and pressure range. Since line-by-line calculations are completely decoupled between 
each wavelength, this reduction corresponds to another three orders of magnitude reduction from line-by-line 
calculations. 

The combined reduced model and cylindrical-slab formulation allow fully coupled flow/radiation calcula- 
tions to be performed for equilibrium air in cylindrical geometries such as that present arc-constrictor flows. 
In addition, the proposed model reduction does not require any additional information in order to “tune” 
the model for a particular set of conditions. This gives the model a truly predictive capability and will allow 
the it to be used in future arc-constrictor design studies. 
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